# Replication Archive for 
# Alexander Coppock, Emily Ekins and David Kirby (2018), 
# "The Long-lasting Effects of Newspaper Op-Eds on Public Opinion",
# Quarterly Journal of Political Science: Vol. 13: No. 1, pp 59-87.

rm(list = ls())

# Uncomment to install
#install.packages(c("tidyverse", "estimatr", "coefplot"))

library(tidyverse)
library(estimatr)
library(coefplot)

# Set working directory to replication archive

load("mturk_opeds_cleaned.rdata")
load("elite_opeds_cleaned.rdata")
load("dv_df.rdata")
source("opeds_source.R")

dvs <-
  c(
    "dv_amtrak_main_w1",
    "dv_climate_main_w1",
    "dv_flat_main_w1",
    "dv_vets_main_w1",
    "dv_wall_main_w1",
    "dv_amtrak_scale_w1",
    "dv_climate_scale_w1",
    "dv_flat_scale_w1",
    "dv_vets_scale_w1",
    "dv_wall_scale_w1"
  )
dv_labels <-
  rep(c("Amtrak", "Climate", "Flat Tax", "Veterans", "Wall Street"), 2)
Zs <-
  rep(c("amtrak", "climate", "paul", "veterans", "wallstreet"), 2)
ages <- levels(mturk_opeds$age_5)

df <- NULL

for (j in 1:length(Zs)) {
  for (k in 1:length(ages)) {
    local_subset <-
      mturk_opeds$Z %in% c("control", Zs[j]) &
      mturk_opeds$age_5 == ages[k]
    
    local_df <- mturk_opeds[local_subset, ]
    local_formula <- formula(paste0(dvs[j], "~Z"))
    
    summary_df <-
      lm_robust(local_formula, data = local_df) %>%
      tidy %>%
      mutate(
        treatment = Zs[j],
        dependent_variable = dv_labels[j],
        dv = dvs[j],
        age = ages[k],
        sample = "Mechanical Turk"
      )
  
    df <- bind_rows(df, summary_df)
  }
}

dvs <- c("dv_amtrak_main_w1", "dv_flat_main_w1", "dv_vets_main_w1", "dv_wall_main_w1",
         "dv_amtrak_scale_w1", "dv_flat_scale_w1", "dv_vets_scale_w1", "dv_wall_scale_w1")
dv_labels <- rep(c("Amtrak", "Flat Tax", "Veterans", "Wall Street"), 2)
Zs <- rep(c("amtrak", "paul", "veterans", "wallstreet"), 2)
ages <- levels(elite_opeds$age_5)


for(j in 1:length(Zs)) {
  for (k in 1:length(ages)) {
    local_subset <-
      elite_opeds$Z %in% c("control", Zs[j]) &
      elite_opeds$age_5 == ages[k]
    local_df <- elite_opeds[local_subset, ]
    local_formula <- formula(paste0(dvs[j], "~Z"))

    summary_df <-
      lm_robust(local_formula, data = local_df) %>%
      tidy %>%
      mutate(
        treatment = Zs[j],
        dependent_variable = dv_labels[j],
        dv = dvs[j],
        age = ages[k],
        sample = "Elites"
      )
    
    df <- bind_rows(df, summary_df)
  }
}


df <-
  df %>%
  filter(term != "(Intercept)") %>%
  mutate(dv_type = ifelse(grepl("scale", outcome), "Scale DV", "Main DV"),
         age = factor(age, levels = c("18 - 29", "30 - 39", "40 - 49", "50 - 59", "60+")))


g <- 
  ggplot(df, aes(y = dependent_variable, x = estimate, group = age, color = age, shape = age)) + 
  geom_point(position = position_dodgev(height = .5), size = 2) +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high),
                 position = position_dodgev(height = .5),
                 height = 0) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
    axis.title = element_blank(),
    legend.key.width = unit(4, "lines"),
    strip.background = element_blank()
  ) +
  ylab("Estimated Treatment Effects") +
  facet_grid(sample ~ dv_type, scales = "free") 

g
